library(foreign)
library(haven)
library(tidyverse)
library(dplyr)
library(plyr)
library(plm)
library(multiwayvcov)
library(lmtest)
library(car)
library(PanelMatch)
library(estimatr)

#############
##FUNCTIONS##
#############

interact.sum <- function(model, vars) {
  coef <- sum(summary(model)$coefficients[vars,1])
  variance <- sum(diag(vcovHC(model, type="HC1")[vars,vars])) +
    2*vcovHC(model, type="HC1")[vars,vars][2,1]
  se <- sqrt(variance)
  output <- data.frame(coef=coef, se=se)
  rownames(output) <- paste(vars, collapse="+")
  return(output)
}

rmse_calculator <- function(i, type) {
  file <- paste0(type,"_model_paper", i, ".rds")
  m <- readRDS(file)[[1]]
  rmse <- sqrt(sum(m$residuals^2)/m$df.residual)
  n <- pdim(m)$nT$N
  return(data.frame(model=i, rmse=rmse, n=n))
}

lags <- function(data, var) {
  n_lags <- length(data[,var])-1
  if(n_lags>0) {
    for(i in 1:n_lags) {
      data[,paste0(var,"_l",i)] <- c(rep(NA,i), data[1:(length(data[,var])-i),var])
    }
  }
  return(data)
}

############
##ANALYSIS##
############

#Set working directory to local output folder#

workingDirectory <- "/Users/ericmcghee/Dropbox/PPIC/UCF/Post-election Covid/academic paper/analysis"

setwd(workingDirectory)

##EVENT STUDIES (Figure 1 & Table A11)##
#All VBM#
m <- plm(pcto ~ vbm_noex + vbm_perm + vbm_app + avr + edr + early + yr2020*cases_mn_centered + pv_marg + active + 
           vbm_all_m3 + vbm_all_m2 + vbm_all_p0 + vbm_all_p1 + vbm_all_p2,
         data=filter(d, reg_all>0), index=c("geoid","year"), model="within", effect="twoways")
keep_vars <- which(!str_detect(rownames(summary(m)$coefficients), "as.factor"))
results <- round(coeftest(m, vcovHC, type="HC1")[keep_vars,], 3)
out8 <- as.data.frame(results) %>% rownames_to_column(var="Variables") %>%
  mutate(index=1:nrow(.)) %>% 
  gather("Type","Value", -c(index, Variables, `t value`, `Pr(>|t|)`)) %>%
  arrange(index)
saveRDS(list(m, out8), "turnout_model_paper8.rds")

out <- data.frame(cycle=-3:2, var="Universal VBM",
                  effect=c(results[str_detect(rownames(results), "vbm_all_m"),1],
                           0,
                           results[str_detect(rownames(results), "vbm_all_p"),1]),
                  effect_se=c(results[str_detect(rownames(results), "vbm_all_m"),2],
                              0,
                              results[str_detect(rownames(results), "vbm_all_p"),2]))

#No-excuse VBM#
m <- plm(pcto ~ vbm_all + vbm_perm + vbm_app + avr + edr + early + yr2020*cases_mn_centered + pv_marg + active + 
           vbm_noex_m3 + vbm_noex_m2 + vbm_noex_p0 + vbm_noex_p1 + vbm_noex_p2,
         data=filter(d, reg_all>0), index=c("geoid","year"), model="within", effect="twoways")
keep_vars <- which(!str_detect(rownames(summary(m)$coefficients), "as.factor"))
results <- round(coeftest(m, vcovHC, type="HC1")[keep_vars,], 3)
out7 <- as.data.frame(results) %>% rownames_to_column(var="Variables") %>%
  mutate(index=1:nrow(.)) %>% 
  gather("Type","Value", -c(index, Variables, `t value`, `Pr(>|t|)`)) %>%
  arrange(index)
saveRDS(list(m, out7), "turnout_model_paper7.rds")

out <- rbind.fill(out, data.frame(cycle=-3:2, var="No-Excuse VBM",
                                  effect=c(results[str_detect(rownames(results), "vbm_noex_m"),1],
                                           0,
                                           results[str_detect(rownames(results), "vbm_noex_p"),1]),
                                  effect_se=c(results[str_detect(rownames(results), "vbm_noex_m"),2],
                                              0,
                                              results[str_detect(rownames(results), "vbm_noex_p"),2])))

#VBM applications#
m <- plm(pcto ~ vbm_all + vbm_noex + vbm_perm + avr + edr + early + yr2020*cases_mn_centered + pv_marg + active + 
           vbm_app_m3 + vbm_app_m2 + vbm_app_p0,
         data=filter(d, reg_all>0), index=c("geoid","year"), model="within", effect="twoways")
keep_vars <- which(!str_detect(rownames(summary(m)$coefficients), "as.factor"))
results <- round(coeftest(m, vcovHC, type="HC1")[keep_vars,], 3)
out6 <- as.data.frame(results) %>% rownames_to_column(var="Variables") %>%
  mutate(index=1:nrow(.)) %>% 
  gather("Type","Value", -c(index, Variables, `t value`, `Pr(>|t|)`)) %>%
  arrange(index)
saveRDS(list(m, out6), "turnout_model_paper6.rds")

out <- rbind.fill(out, data.frame(cycle=-3:0, var="VBM Application",
                                  effect=c(results[str_detect(rownames(results), "vbm_app_m"),1],
                                           0,
                                           results[str_detect(rownames(results), "vbm_app_p"),1]),
                                  effect_se=c(results[str_detect(rownames(results), "vbm_app_m"),2],
                                              0,
                                              results[str_detect(rownames(results), "vbm_app_p"),2])))

out %>% ggplot(aes(cycle, effect)) + 
  scale_x_continuous(breaks=seq(-3,2)) +
  scale_y_continuous(breaks=seq(-0.08,0.08,by=0.02), labels=scales::label_percent()) +
  coord_cartesian(ylim=c(-0.08,0.08)) +
  labs(x="Cycles until/since (-/+) reform", y="Turnout deviation (Baseline t-1)") +
  geom_hline(yintercept=0, color="gray55") +
  geom_vline(xintercept=-0.5, linetype="dashed") +
  geom_line() +
  geom_point() + 
  geom_errorbar(aes(ymin=effect-2*effect_se, ymax=effect+2*effect_se), width=0.1) +
  facet_wrap(~var)

ggsave("figure_1.png")

##DIFFERENCE-DIFFERENCES (Figure 3 & Table A8)##
#Table A8, Model 5: diff-in-diff w/ 2020 interaction#
m <- plm(pcto ~ yr2020*(vbm_all + vbm_noex + vbm_app), 
         filter(d, reg_all>0), index=c("geoid","year"), model="within", effect="twoways")
keep_vars <- which(!str_detect(rownames(summary(m)$coefficients), "as.factor"))
out5 <- as.data.frame(round(coeftest(m, vcovHC, type="HC1")[keep_vars,], 3)) %>% rownames_to_column(var="Variables") %>%
  mutate(index=1:nrow(.)) %>% 
  gather("Type","Value", -c(index, Variables, `t value`, `Pr(>|t|)`)) %>%
  arrange(index)
saveRDS(list(m, out5), "turnout_model_paper5.rds")

#Table A8, Model 4: diff-in-diff w/ 2020 interaction, county trends, & other reforms#
m <- plm(pcto ~ yr2020*(vbm_all + vbm_noex + vbm_app + avr + cases_mn_centered + pv_marg) + vbm_perm + edr + early + active, 
         filter(d, reg_all>0), index=c("geoid","year"), model="within", effect="twoways")
keep_vars <- which(!str_detect(rownames(summary(m)$coefficients), "as.factor"))
out4 <- as.data.frame(round(coeftest(m, vcovHC, type="HC1")[keep_vars,], 3)) %>% rownames_to_column(var="Variables") %>%
  mutate(index=1:nrow(.)) %>% 
  gather("Type","Value", -c(index, Variables, `t value`, `Pr(>|t|)`)) %>%
  arrange(index)
saveRDS(list(m, out4), "turnout_model_paper4.rds")

#Table A8, Model 3: diff-in-diff w/ 2020 interaction, county trends, & shorter time frame#
m <- plm(pcto ~ yr2020*(vbm_all + vbm_noex + vbm_app + avr + cases_mn_centered + pv_marg) + vbm_perm + edr + early + active +
           (as.integer(as.factor(year))-1)*as.factor(geoid), 
         data=filter(d, year>=2008, reg_all>0), index=c("geoid","year"), model="within", effect="twoways")
keep_vars <- which(!str_detect(rownames(summary(m)$coefficients), "as.factor"))
out3 <- as.data.frame(round(coeftest(m, vcovHC, type="HC1")[keep_vars,], 3)) %>% rownames_to_column(var="Variables") %>%
  mutate(index=1:nrow(.)) %>% 
  gather("Type","Value", -c(index, Variables, `t value`, `Pr(>|t|)`)) %>%
  arrange(index)
saveRDS(list(m, out3), "turnout_model_paper3.rds")

#Table A8, Model 2: diff-in-diff w/ 2020 interaction, county trends, & other reforms#
m <- plm(pcto ~ yr2020*(vbm_all + vbm_noex + vbm_app + avr + cases_mn_centered + pv_marg) + vbm_perm + edr + early + active + 
           yrcount*as.factor(geoid), 
         filter(d, reg_all>0), index=c("geoid","year"), model="within", effect="twoways")
keep_vars <- which(!str_detect(rownames(summary(m)$coefficients), "as.factor"))
out2 <- as.data.frame(round(coeftest(m, vcovHC, type="HC1")[keep_vars,], 3)) %>% rownames_to_column(var="Variables") %>%
  mutate(index=1:nrow(.)) %>% 
  gather("Type","Value", -c(index, Variables, `t value`, `Pr(>|t|)`)) %>%
  arrange(index)
saveRDS(list(m, out2), "turnout_model_paper2.rds")

#Table A8, Model 1: diff-in-diff w/ 2020 interaction, county trends, & percent vbm interaction#
m <- plm(pcto ~ yr2020*(vbm_all + pcvbm + vbm_noex + vbm_app + avr + cases_mn_centered + pv_marg) + vbm_perm + edr + early + active + 
           yrcount*as.factor(geoid), 
         filter(d, reg_all>0), index=c("geoid","year"), model="within", effect="twoways")
keep_vars <- which(!str_detect(rownames(summary(m)$coefficients), "as.factor"))
out1 <- as.data.frame(round(coeftest(m, vcovHC, type="HC1")[keep_vars,], 3)) %>% rownames_to_column(var="Variables") %>%
  mutate(index=1:nrow(.)) %>% 
  gather("Type","Value", -c(index, Variables, `t value`, `Pr(>|t|)`)) %>%
  arrange(index)
saveRDS(list(m, out1), "turnout_model_paper1.rds")

out1 <- readRDS("turnout_model_paper1.rds")[[2]]
out2 <- readRDS("turnout_model_paper2.rds")[[2]]
out3 <- readRDS("turnout_model_paper3.rds")[[2]]
out4 <- readRDS("turnout_model_paper4.rds")[[2]]
out5 <- readRDS("turnout_model_paper5.rds")[[2]]
out6 <- readRDS("turnout_model_paper6.rds")[[2]]
out7 <- readRDS("turnout_model_paper7.rds")[[2]]
out8 <- readRDS("turnout_model_paper8.rds")[[2]]

out.all <- full_join(out1, out2, by=c("Variables","Type"), suffix=c(".1",".2")) %>%
  full_join(out3, by=c("Variables","Type"), suffix=c("",".3")) %>%
  full_join(out4, by=c("Variables","Type"), suffix=c("",".4")) %>%
  full_join(out5, by=c("Variables","Type"), suffix=c("",".5")) %>%
  full_join(out6, by=c("Variables","Type"), suffix=c("",".6")) %>%
  full_join(out7, by=c("Variables","Type"), suffix=c("",".7")) %>%
  full_join(out8, by=c("Variables","Type"), suffix=c("",".8"))

out.all <- out.all[,!str_detect(names(out.all), "t value|Pr|index|Type")]

write.csv(out.all, "turnout_results_table.csv")

#data for figure 3#
figure.results <- data.frame(var=c("Universal VBM","Universal VBM among Precinct Voters", "No-Excuse VBM", "VBM Application"), 
                             est=c(out.all$Value.2[1], out.all$Value.1[c(1,5,7)]),
                             se=c(out.all$Value.2[2], out.all$Value.1[c(2,6,8)]))

m <- readRDS("turnout_model_paper2.rds")[[1]]
vars <- c("vbm_all","yr2020TRUE:vbm_all")
int <- round(interact.sum(m, vars), 3)
figure.results <- rbind.fill(figure.results,
                             data.frame(var="Universal VBM in 2020",
                                        est=int[[1]],
                                        se=int[[2]]))

m <- readRDS("turnout_model_paper1.rds")[[1]]
vars <- c("vbm_all","yr2020TRUE:vbm_all")
int <- round(interact.sum(m, vars), 3)
figure.results <- rbind.fill(figure.results,
                             data.frame(var="Universal VBM among Precinct Voters in 2020",
                                        est=int[[1]],
                                        se=int[[2]]))
vars <- c("vbm_noex","yr2020TRUE:vbm_noex")
int <- round(interact.sum(m, vars), 3)
figure.results <- rbind.fill(figure.results,
                             data.frame(var="No-Excuse VBM in 2020",
                                        est=int[[1]],
                                        se=int[[2]]))
vars <- c("vbm_app","yr2020TRUE:vbm_app")
int <- round(interact.sum(m, vars), 3)
figure.results <- rbind.fill(figure.results,
                             data.frame(var="VBM Application in 2020",
                                        est=int[[1]],
                                        se=int[[2]]))
saveRDS(figure.results, "data_for_figure_3.rds")

#figure 3#
f <- readRDS("data_for_figure_3.rds") %>%
  mutate(yr2020=c(rep("Before 2020",4),rep("In 2020",4)),
         var=c(var[1:4],var[1:4]),
         index=4:1) %>%
  mutate(var=factor(var, levels=unique(var)[order(unique(index))]))

f %>% ggplot(aes(est, var)) + 
  scale_x_continuous(breaks=seq(-0.1,0.15,by=0.05), labels=scales::label_percent()) +
  coord_cartesian(xlim=c(-0.1,0.15)) +
  labs(y="", x="Change in total turnout") +
  geom_vline(xintercept=0, linetype="dashed") +
  geom_point() + 
  geom_errorbar(aes(xmin=est-2*se, xmax=est+2*se), width=0.1) +
  facet_grid(yr2020~.)

ggsave("figure_3.png")

#RMSEs for tables#
model_rmses <- ldply(lapply(1:8, function(i,t) rmse_calculator(i,t), "turnout"))
saveRDS(model_rmses, "turnout_model_fits.rds")

##LAGGED DV MODEL (Tables 1 & A4)##
d2 <- arrange(d, stpost, geoid, year) %>%
  filter(!is.na(geoid)) %>%
  ddply(.(geoid), function(x,y) lags(x,y), "pcto")

m <- lm_robust(pcto ~ vbm_all + vbm_noex + vbm_app + avr + cases_mn_centered + pv_marg + active +
         pcto_l5 + pcto_l4 + pcto_l3 + pcto_l2 + pcto_l1,
        data=filter(d2, year==2020), clusters=stpost)
summary(m)

out.ldv <- as.data.frame(round(coeftest(m)[1:nrow(coeftest(m)),], 3)) %>% rownames_to_column(var="Variables") %>%
  mutate(index=1:nrow(.)) %>% 
  gather("Type","Value", -c(index, Variables, `t value`, `Pr(>|t|)`)) %>%
  arrange(index)
saveRDS(list(m, out.ldv), "turnout_ldv_paper.rds")

rmse <- sqrt(sum(m$residuals^2)/m$df.residual)
n <- summary(m)$df[1] + summary(m)$df[2]
saveRDS(data.frame(model="ldv", rmse=rmse, n=n), "turnout_ldv_fit.rds")

##MATCHING MODELS (Tables 1 & A3)##
#Universal VBM#
d <- mutate(d, geoid2=as.integer(geoid), yrcount2=as.integer(yrcount+1)) %>%
  ddply(.(geoid), transform, pre_2020=any(vbm_all[year<2020]))

DisplayTreatment(unit.id = "geoid2",
                 time.id = "yrcount2", legend.position = "none",
                 xlab = "Year", ylab = "County Code",
                 treatment = "vbm_all", data = filter(d, !(pre_2020)), dense.plot=T,
                 hide.y.axis.label = T)

pm.out <- PanelMatch(lag = 5, time.id = "yrcount2", unit.id = "geoid2", 
                     treatment = "vbm_all", refinement.method = "none", # no refinement 
                     data = filter(d, !(pre_2020)), 
                     match.missing = TRUE, 
                     covs.formula = ~ I(lag(pcto, 1:5)) +  I(lag(pv_marg, 1:5)) +
                       vbm_noex + vbm_app + avr + cases_mn_centered + pv_marg, 
                     size.match = 5, qoi = "att" , outcome.var = "pcto",
                     lead = 0, forbid.treatment.reversal = FALSE, 
                     use.diagonal.variance.matrix = TRUE)

pm.balance.vbm_all.null <- get_covariate_balance(pm.out$att,
                                                  data = filter(d, !(pre_2020)),
                                                  covariates = c("cases_mn_centered","pv_marg","pcto"),
                                                  plot = FALSE)
pm.balance.vbm_all.null
pm.balance.vbm_all.null["t_0","pcto"] <- NA
pm.balance <- cbind(data.frame(dv="vbm_all", model="null"), 
                    value=apply(pm.balance.vbm_all.null, 2, function(x) mean(x, na.rm=T)))

pm.out <- PanelMatch(lag = 5, time.id = "yrcount2", unit.id = "geoid2", 
                     treatment = "vbm_all", refinement.method = "mahalanobis", # use Mahalanobis distance 
                     data = filter(d, !(pre_2020)), 
                     match.missing = TRUE, 
                     covs.formula = ~ I(lag(pcto, 1:5)) +  I(lag(pv_marg, 1:5)) +
                       vbm_noex + vbm_app + avr + cases_mn_centered + pv_marg, 
                     size.match = 5, qoi = "att" , outcome.var = "pcto",
                     lead = 0, forbid.treatment.reversal = FALSE, 
                     use.diagonal.variance.matrix = TRUE)

pm.balance.vbm_all.match <- get_covariate_balance(pm.out$att,
                                                  data = filter(d, !(pre_2020)),
                                                  covariates = c("cases_mn_centered","pv_marg","pcto"),
                                                  plot = FALSE)
pm.balance.vbm_all.match
pm.balance.vbm_all.match["t_0","pcto"] <- NA
pm.balance <- rbind.fill(pm.balance, cbind(data.frame(dv="vbm_all", model="matched"), 
                                           value=apply(pm.balance.vbm_all.match, 2, function(x) mean(x, na.rm=T))))

pe.out.vbm_all.match <- PanelEstimate(sets=pm.out, data=filter(d, !(pre_2020)))
summary(pe.out.vbm_all.match)

#No-excuse VBM#
d <- mutate(d, geoid2=as.integer(geoid), yrcount2=as.integer(yrcount+1)) %>%
  ddply(.(geoid), transform, pre_2020=any(vbm_noex[year<2020]))

pm.out <- PanelMatch(lag = 5, time.id = "yrcount2", unit.id = "geoid2", 
                     treatment = "vbm_noex", refinement.method = "none", # no refinement 
                     data = filter(d, !(pre_2020)), 
                     match.missing = TRUE, 
                     covs.formula = ~ I(lag(pcto, 1:5)) +  I(lag(pv_marg, 1:5)) +
                       vbm_all + vbm_app + avr + cases_mn_centered, 
                     size.match = 5, qoi = "att" , outcome.var = "pcto",
                     lead = 0, forbid.treatment.reversal = FALSE, 
                     use.diagonal.variance.matrix = TRUE)

pm.balance.vbm_noex.null <- get_covariate_balance(pm.out$att,
                                                  data = filter(d, !(pre_2020)),
                                                  covariates = c("cases_mn_centered","pv_marg","pcto"),
                                                  plot = FALSE)
pm.balance.vbm_noex.null
pm.balance.vbm_noex.null["t_0","pcto"] <- NA
pm.balance <- rbind.fill(pm.balance, cbind(data.frame(dv="vbm_noex", model="null"), 
                                           value=apply(pm.balance.vbm_noex.null, 2, function(x) mean(x, na.rm=T))))

pm.out <- PanelMatch(lag = 5, time.id = "yrcount2", unit.id = "geoid2", 
                     treatment = "vbm_noex", refinement.method = "mahalanobis", # use Mahalanobis distance 
                     data = filter(d, !(pre_2020)), 
                     match.missing = TRUE, 
                     covs.formula = ~ I(lag(pcto, 1:5)) +  I(lag(pv_marg, 1:5)) +
                       vbm_all + vbm_app + avr + cases_mn_centered, 
                     size.match = 5, qoi = "att" , outcome.var = "pcto",
                     lead = 0, forbid.treatment.reversal = FALSE, 
                     use.diagonal.variance.matrix = TRUE)

pm.balance.vbm_noex.match <- get_covariate_balance(pm.out$att,
                                                   data = filter(d, !(pre_2020)),
                                                   covariates = c("cases_mn_centered","pv_marg","pcto"),
                                                   plot = FALSE)
pm.balance.vbm_noex.match
pm.balance.vbm_noex.match["t_0","pcto"] <- NA
pm.balance <- rbind.fill(pm.balance, cbind(data.frame(dv="vbm_noex", model="matched"), 
                                           value=apply(pm.balance.vbm_noex.match, 2, function(x) mean(x, na.rm=T))))

pe.out.vbm_noex.match <- PanelEstimate(sets=pm.out, data=filter(d, !(pre_2020)))
summary(pe.out.vbm_noex.match)

#VBM application#
d <- mutate(d, geoid2=as.integer(geoid), yrcount2=as.integer(yrcount+1)) %>%
  ddply(.(geoid), transform, pre_2020=any(vbm_app[year<2020]))

pm.out <- PanelMatch(lag = 5, time.id = "yrcount2", unit.id = "geoid2", 
                     treatment = "vbm_app", refinement.method = "none", # no refinement 
                     data = filter(d, !(pre_2020)), 
                     match.missing = TRUE, 
                     covs.formula = ~ I(lag(pcto, 1:5)) +  I(lag(pv_marg, 1:5)) +
                       vbm_all + vbm_noex + avr + cases_mn_centered, 
                     size.match = 5, qoi = "att" , outcome.var = "pcto",
                     lead = 0, forbid.treatment.reversal = FALSE, 
                     use.diagonal.variance.matrix = TRUE)

pm.balance.vbm_app.null <- get_covariate_balance(pm.out$att,
                                                 data = filter(d, !(pre_2020)),
                                                 covariates = c("cases_mn_centered","pv_marg","pcto"),
                                                 plot = FALSE)
pm.balance.vbm_app.null
pm.balance.vbm_app.null["t_0","pcto"] <- NA
pm.balance <- rbind.fill(pm.balance, cbind(data.frame(dv="vbm_app", model="null"), 
                                           value=apply(pm.balance.vbm_app.null, 2, function(x) mean(x, na.rm=T))))

pm.out <- PanelMatch(lag = 5, time.id = "yrcount2", unit.id = "geoid2", 
                     treatment = "vbm_app", refinement.method = "mahalanobis", # use Mahalanobis distance 
                     data = filter(d, !(pre_2020)), 
                     match.missing = TRUE, 
                     covs.formula = ~ I(lag(pcto, 1:5)) +  I(lag(pv_marg, 1:5)) +
                       vbm_all + vbm_noex + avr + cases_mn_centered, 
                     size.match = 5, qoi = "att" , outcome.var = "pcto",
                     lead = 0, forbid.treatment.reversal = FALSE, 
                     use.diagonal.variance.matrix = TRUE)

pm.balance.vbm_app.match <- get_covariate_balance(pm.out$att,
                                                  data = filter(d, !(pre_2020)),
                                                  covariates = c("cases_mn_centered","pv_marg","pcto"),
                                                  plot = FALSE)
pm.balance.vbm_app.match
pm.balance.vbm_app.match["t_0","pcto"] <- NA
pm.balance <- rbind.fill(pm.balance, cbind(data.frame(dv="vbm_app", model="match"), 
                                           value=apply(pm.balance.vbm_app.match, 2, function(x) mean(x, na.rm=T))))

pe.out.vbm_app.match <- PanelEstimate(sets=pm.out, data=filter(d, !(pre_2020)))
summary(pe.out.vbm_app.match)

pm.balance <- rbind.fill(cbind(data.frame(var=colnames(pm.balance.vbm_all.null), dv="vbm_all", model="null"), 
                               value=apply(pm.balance.vbm_all.null, 2, function(x) mean(x, na.rm=T))),
                         cbind(data.frame(var=colnames(pm.balance.vbm_all.match), dv="vbm_all", model="matched"), 
                               value=apply(pm.balance.vbm_all.match, 2, function(x) mean(x, na.rm=T))),
                         cbind(data.frame(var=colnames(pm.balance.vbm_noex.null), dv="vbm_noex", model="null"), 
                               value=apply(pm.balance.vbm_noex.null, 2, function(x) mean(x, na.rm=T))),
                         cbind(data.frame(var=colnames(pm.balance.vbm_noex.match), dv="vbm_noex", model="matched"), 
                               value=apply(pm.balance.vbm_noex.match, 2, function(x) mean(x, na.rm=T))),
                         cbind(data.frame(var=colnames(pm.balance.vbm_app.null), dv="vbm_app", model="null"), 
                               value=apply(pm.balance.vbm_app.null, 2, function(x) mean(x, na.rm=T))),
                         cbind(data.frame(var=colnames(pm.balance.vbm_app.match), dv="vbm_app", model="matched"), 
                               value=apply(pm.balance.vbm_app.match, 2, function(x) mean(x, na.rm=T))))

write.csv(spread(pm.balance, model, value) %>% arrange(dv, var), "turnout_matching_balance.csv")

pe.out <- rbind.fill(data.frame(summary(pe.out.vbm_all.match)$summary),
                     data.frame(summary(pe.out.vbm_noex.match)$summary),
                     data.frame(summary(pe.out.vbm_app.match)$summary))

write.csv(pe.out, "turnout_matching_estimates.csv")


